#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <opencv2/opencv.hpp>
//#include <Eigen/Eigen>
#include "fstream"
#include "dragon.h"
//using namespace std;
//using namespace cv;
//using namespace Eigen;
int width = 32, height = 32;
const int samplesize = 50000;
void initCubemapPoint(cv::Mat &cubeface)
{
    //0为正z，1为-z，2为正x，3为负x，4为正y（x朝上），5为负y,都是从里面看，从上到下，从左到右
    int id=0;
    //+z
    for(int i=0;i<height;i++){
        for(int j=0;j<width;j++){
            cubeface.at<float>(id, 0) = 16-j;
            cubeface.at<float>(id, 1) = 16-i;
            cubeface.at<float>(id, 2) = 16;
            id++;
        }
    }
    //-z
    for(int i=0;i<height;i++){
        for(int j=0;j<width;j++){
            cubeface.at<float>(id, 0) = -16+j;
            cubeface.at<float>(id, 1) = 16-i;
            cubeface.at<float>(id, 2) = -16;
            id++;
        }
    }
    //+x
    for(int i=0;i<height;i++){
        for(int j=0;j<width;j++){
            cubeface.at<float>(id, 0) = 16;
            cubeface.at<float>(id, 1) = 16-i;
            cubeface.at<float>(id, 2) = -16+j;
            id++;
        }
    }
    //-x
    for(int i=0;i<height;i++){
        for(int j=0;j<width;j++){
            cubeface.at<float>(id, 0) = -16;
            cubeface.at<float>(id, 1) = 16-i;
            cubeface.at<float>(id, 2) = 16-j;
            id++;
        }
    }
    //+y
    for(int i=0;i<height;i++){
        for(int j=0;j<width;j++){
            cubeface.at<float>(id, 0) = 16-i;
            cubeface.at<float>(id, 1) = 16;
            cubeface.at<float>(id, 2) = 16-j;
            id++;
        }
    }
    //-y
    for(int i=0;i<height;i++){
        for(int j=0;j<width;j++){
            cubeface.at<float>(id, 0) = 16-i;
            cubeface.at<float>(id, 1) = -16;
            cubeface.at<float>(id, 2) = -16+j;
            id++;
        }
    }

}
void calculateDot(cv::Mat& dot, const cv::Mat& cubeface)
{
    for(int i=0;i<dot.rows;i++){
        cv::Mat curpoint(1,3,CV_32FC1, normals[i]);//当前点
        //std::cout<<curpoint<<std::endl;
        normalize(curpoint,curpoint);
        //std::cout<<curpoint<<std::endl;
        cv::Mat point_mat = repeat(curpoint, 6*32*32,1);
        //std::cout<<point_mat<<std::endl;
        //光方向
        cv::Mat lightdirect(6*32*32,3,CV_32FC1);
        lightdirect = cubeface.clone();
        //std::cout<<lightdirect<<std::endl;
        for(int j=0; j<lightdirect.rows; j++){
            normalize(lightdirect.row(j), lightdirect.row(j));
        }
        point_mat = point_mat.mul(lightdirect);
        //std::cout<<point_mat<<std::endl;
        transpose(point_mat.col(0)+point_mat.col(1)+point_mat.col(2),dot.row(i));
        //std::cout<<dot.row(i)<<std::endl;cur
        if(i%5000==0)
            std::cout<<"dot"<<i<<std::endl;
    }
}
void calFresnelTrans(cv::Mat& fresnel, const cv::Mat& dot)
{
    //下面计算一脸懵逼的fresnel transmittance
    //wiki说transmittance = n2*cos(theta_t)/(n1*cos(theta_i)
    //yita = n2/n1=sin(theta_i)/sin(theta_t), n2代表物体材质，n1是空气，theta_t是折射角， theta_i是入射角
    float yita = 1.3;
    cv::Mat theta_t(dot.rows, dot.cols, CV_32FC1);
    theta_t = cv::abs(dot)/yita;//因为法向量和入射光都是单位向量，所以点乘结果的绝对值就是入射角的sin，此步求出折射角的sin
    theta_t = 1- theta_t.mul(theta_t);//求出折射角的cos
    for(int i=0; i<theta_t.rows;i++){
        for(int j=0;j<theta_t.cols;j++){
            theta_t.at<float>(i,j) = std::sqrt(theta_t.at<float>(i,j));
        }
    }
    cv::Mat cur = dot.mul(dot);
    for(int i=0; i<cur.rows;i++){
        for(int j=0;j<cur.cols;j++){
            cur.at<float>(i,j) = std::sqrt(1-cur.at<float>(i,j));
        }
    }
    cv::divide((theta_t.mul(cur)*4).mul(cur), (cur+yita*theta_t).mul(cur+yita*theta_t), theta_t);
    cv::divide(yita*theta_t, cur,fresnel);//amplitude ratio 什么鬼

    std::cout<<"hh"<<std::endl;
}
const float PI = 3.1415;
float getR(const float& r, const float& z_r, const float& A, const float& alpha2,const float& sigma_tr)
{
    float d_r = sqrt(r*r + z_r*z_r),z_v = z_r*(1+4*A/3), d_v = sqrt(r*r +z_v*z_v);
    float R;
    R = alpha2/(4*PI)*(z_r*(sigma_tr+1/d_r)*exp(-sigma_tr*d_r)/d_r/d_r +
                       z_v*(sigma_tr+1/d_v)*exp(-sigma_tr*d_v)/d_v/d_v);
    return R;
}
void getAllR(cv::Mat& R)
{
    float yita = 1.5;
    int samplesize = 50000;
    float sigma_s = 0.75, sigma_a = 0.02;
    float g=0.25;//散射角平均cos
    float sigma_s2 = (1-g)*sigma_s, sigma_t2 = sigma_a+sigma_s2,
            alpha2 = sigma_s2/sigma_t2, sigma_tr = sqrt(3*sigma_a*sigma_t2), z_r = 1/sigma_t2,
            F_dr = -1.440/yita*yita+0.710/yita+0.668+0.0636*yita, A = (1+F_dr)*(1-F_dr);
    for(int i=0;i<samplesize;i++){
        for(int j=0;j<samplesize;j++){
            float cur_r = sqrt(pow(vertices[i][0]-vertices[j][0],2)+
                    pow(vertices[i][1]-vertices[j][1],2)+
                    pow(vertices[i][2]-vertices[j][2],2));
            R.at<float>(i,j) = getR(cur_r,z_r,A,alpha2,sigma_tr);
        }
    }
}
cv::Mat Wave(cv::Mat& img)
{
     int Height = img.cols;
     int Width = img.rows;
     int depth = 3; //定义分解深度
     int depthcount = 1;
     cv::Mat tmp = cv::Mat::ones(Width, Height, CV_32FC1);
     cv::Mat wavelet = cv::Mat::ones(Width, Height, CV_32FC1);
     cv::Mat imgtmp = img.clone();
     imgtmp.convertTo(imgtmp, CV_32FC1);
     while (depthcount<=depth){
         Width = img.rows / depthcount;
         Height = img.cols / depthcount;
         for (int i = 0; i < Width; i++){
             for (int j = 0; j < Height / 2; j++){
                 tmp.at<float>(i, j) = (imgtmp.at<float>(i, 2 * j) + imgtmp.at<float>(i, 2 * j + 1)) / 2;
                 tmp.at<float>(i, j + Height / 2) = (imgtmp.at<float>(i, 2 * j) - imgtmp.at<float>(i, 2 * j + 1)) / 2;
             }
         }
         for (int i = 0; i < Width / 2; i++){
             for (int j = 0; j < Height; j++){
                 wavelet.at<float>(i, j) = (tmp.at<float>(2 * i, j) + tmp.at<float>(2 * i + 1, j)) / 2;
                 wavelet.at<float>(i + Width / 2, j) = (tmp.at<float>(2 * i, j) - tmp.at<float>(2 * i + 1, j)) / 2;
             }
         }
         imgtmp = wavelet;
         depthcount++;
     }
     wavelet.convertTo(wavelet, CV_8UC1);
    // wavelet += 50; //图像暗度过低，所以这里我加了50
     imshow("jpg", wavelet);
     return wavelet;
}
void WaveE(cv::Mat& E)
{
    for(int i=0;i<E.rows;i++){
        for(int j=0;j<6;j++){
            cv::Mat cur = E.row(i).colRange(32*32*j,32*32*(j+1));
            cur.reshape(0, 32);
            Wave(cur).copyTo(cur);
            cur.reshape(0,1);
            cur.copyTo(E.row(i).colRange(32*32*j, 32*32*(j+1)));
        }
    }

}
std::vector<std::string> filename;
void InitVisibility(cv::Mat& visibility)
{
    filename.push_back("./v/visibility1.txt");
    filename.push_back("./v/visibility2.txt");
    filename.push_back("./v/visibility3.txt");
    filename.push_back("./v/visibility4.txt");
    filename.push_back("./v/visibility5.txt");
    filename.push_back("./v/visibility6.txt");
    filename.push_back("./v/visibility7.txt");
    filename.push_back("./v/visibility8.txt");
    filename.push_back("./v/visibility9.txt");
    filename.push_back("./v/visibility10.txt");
    int row = visibility.rows, col = visibility.cols;
    for(int k=0;k<10;k++){
        std::string name = filename[k];
        std::ifstream in(name.c_str());
        for(int i=0;i<5000;i++){
            uchar* p = visibility.ptr<uchar>(k*5000+i);
            for(int j=0;j<visibility.cols;j++){
                int cur;
                in>>cur;
                *p = cur;
                p++;
            }
        }
    }
}

int main (int argc,  char *argv[])
{
//    cv::Mat src0 = imread("back.png",0);
//    cv::Mat src(src0.rows, src0.cols,CV_32FC1);
//    src0.copyTo(src);
//    Wave(src);


    //std::cout<<visibility.row(samplesize-1)<<std::endl;
    const int cubesize = 32;
    cv::Mat cubeface(cubesize*cubesize*6 ,3, CV_32FC1);
    cv::Mat dot(samplesize, cubesize*cubesize*6, CV_32FC1);

    /*initCubemapPoint(cubeface);//计算光源坐标
    std::cout<<"cubeface is finished"<<std::endl;
    calculateDot(dot, cubeface);
    std::cout<<"calculateDot is finished"<<std::endl;
    std::ofstream out_dot("dot.txt", std::ios::binary);
    out_dot.write((char*)dot.data, dot.rows * dot.cols * sizeof(float));
    std::cout<<"output dot finish"<<std::endl;*/
//    std::cout<<"open before"<<std::endl;

    std::ifstream in("dot.txt", std::ios::binary);
    std::cout<<"open"<<std::endl;

    in.read((char*)dot.data, dot.rows * dot.cols * sizeof(float));
        cubeface.release();
    cv::Mat fresnel(dot.rows, dot.cols, CV_32FC1);
    calFresnelTrans(fresnel, dot);
    std::cout<<"calFresnelTrans is finished"<<std::endl;

    cv::Mat visibility(samplesize, 6*32*32, CV_8UC1);
    InitVisibility(visibility);

    cv::Mat E = cv::Mat::zeros(fresnel.rows, fresnel.cols, CV_32FC1);
    for(int i=0;i<E.rows;i++){
        float* p = E.ptr<float>(i);
        uchar* p_v = visibility.ptr<uchar>(i);
        float* p_f = fresnel.ptr<float>(i);
        float* p_d = dot.ptr<float>(i);
        for(int j=0;j<E.cols;j++){
            if(!(*p_v))
                *p = (*p_f)*(*p_d);
            p_v++;
            p++;
            p_f++;
            p_d++;
        }
    }
    std::cout<<"calE is finished"<<std::endl;
    std::ofstream out("E.txt",std::ios::binary);
    out.write((char*)E.data, E.rows * E.cols * sizeof(float));
    std::cout<<"output to file is finished"<<std::endl;
    /*WaveE(E);//对E做小波变换
    E.colRange(0,128).copyTo(E);
    cv::Mat R(samplesize,samplesize,CV_32FC1);
    getAllR(R);*/
}

